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ABSTRACT 

Instabilities associated with fluid handling and operation in liquid rocket propulsion systems and 
test facilities usually manifest themselves as structural vibrations or some form of structural damage. 
While the source of the instability is directly related to the performance of a component such as a 
turbopump, valve or a flow control element, the associated pressure fluctuations as they propagate 
through the system have the potential to amplify and resonate with natural modes of the structural 
elements and components of the system. In this paper, the authors have developed an innovative multi- 
level approach that involves analysis at the component and systems level. The primary source of the 
unsteadiness is modeled with a high-fidelity hybrid RANS/LES based CFD methodology that has been 
previously used to study instabilities in feed systems. This high fidelity approach is used to quantify the 
instability and understand the physics associated with the instability. System response to the driving 
instability is determined through a transfer matrix approach wherein the incoming and outgoing pressure 
and velocity fluctuations are related through a transfer (or transmission) matrix. The coefficients of the 
transfer matrix for each component (i.e. valve, pipe, orifice etc.) are individually derived from the flow 
physics associated with the component. A demonstration case representing a test loop/test facility 
comprised of a network of elements is constructed with the transfer matrix approach and the amplification 
of modes analyzed as the instability propagates through the test loop. 

INTRODUCTION 

Rocket engine and component testing facilities are usually comprised of an assembly of elements 
such as cavitating venturis, orifices, turbopumps, turning ducts and control valves that are networked 
together in a configuration that meets design goals and testing requirements. The complexity in the 
flowpath of such configurations, however, has the potential of introducing flow induced instabilities due to 
vortex shedding, cavitation, turbulence and large scale flow separation. Such flow phenomena are 
usually accompanied by pressure fluctuations that can interact with test system components as well as 
test articles such as pumps and combustors. Alternately, the pressure fluctuations can also arise from 
the transient stages of testing such as during engine start-up, shutdown or during sudden changes in 
thrust levels. The pressure fluctuations, whether due to flow instability or transient testing, may excite 
certain structural modes and can lead to system wide resonance, thereby compromising the safety of the 
test stands and leading to premature shutdown of tests. Such deleterious structural vibrations coupled 
with large scale pressure fluctuations have been known to lead to system malfunction as has been seen 
in the case of the Delta 4 system. A good example of this coupling is represented by a commonly known 
phenomenon called the “POGO” instability where a closed loop interaction develops between longitudinal 
structural modes of vibration and the instability dynamics of the propulsion system (feed lines and tanks). 
The instability propagates into the pumps where fluctuating inlet pressures result in cavitation leading to 
oscillating thrust, which in turn amplifies the longitudinal structural modes that set the instability in motion 
in the first place. 
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From an analysis perspective, the instability mechanisms seen in testing facilities can be 
classified into three distinct classes: (a) acoustic/hydrodynamic/fluid dynamic instabilities [1] that are 
predominantly attributed to fundamental flow physical mechanisms such as vortex shedding, turbulence, 
etc. largely associated with complex structural configurations such as valve housing, plug shapes, 
manifolds and bends in the piping system, (b) flow transients associated with valve timing and valve 
scheduling [2], Such transients can play an important role especially during startup\shutdown and valve 
response can be critical for safe and reliable operation. The water hammer instability, being a prime 
example, of rapid shutdown of a valve leading to large amplitude pressure fluctuations in the system, (c) 
multi-phase instabilities such as cavitation related instability mechanisms [3] that become especially 
important with cryogenic working fluids due to reduced liquid to vapor density ratios and strong 
evaporative cooling effects resulting in local temperature fluctuations that can couple with the primary 
cavitation instability. 

It has become commonplace to use specialized CFD based numerical tools to identify the above- 
mentioned instability modes and the dominant frequencies associated with them. For example, high 
fidelity computational analyses have been performed for a series of structurally complex valve systems 
such as the pressure regulator valve, split-body valve and globe valves. Such simulations have been 
usually carried out with an advanced hybrid RANS-LES turbulence framework [4] that has permitted 
accurate representation of a rich spectrum of spatial and temporal scales very efficiently. Moreover, a 
novel valve motion computational technique coupled with a grid movement algorithm has been used to 
simulate transient motion of valves during operation. This has permitted analyses of valve scheduling 
procedures and prediction of transients during valve operations. Coupled fluid-structure simulations of 
valve systems have also been undertaken to predict structural response to the governing flow instability 
[5]. In such simulations the effects of the valve packing/dampers have also been included and 
consequently its feedback and the fluid response have been captured. 

For application to flows with multi-phase instabilities the CFD framework was enhanced with sub- 
models for unsteady cavitation and real fluid thermodynamic property routines [6]. These developments 
have facilitated the simulation of two-phase regimes in cryogenic fluids where phase change and 
temperature fluctuations are strongly coupled to pressure and flow variations. Analysis of cavitation 
instabilities such as traveling cloud cavitation have been performed for control elements in test loops such 
as venturis and orifices as well as rotating cavitation for high pressure inducers/turbopumps. 

While the high fidelity CFD simulations discussed above have been successful at identifying the 
physical mechanisms responsible for the fluid instabilities and the fundamental modes associated with 
these instabilities, uncertainties persist on the coupling of the instability modes with the functioning of 
other elements in the test loop, their response to the instabilities, and consequently the resulting structural 
excitation and vibrations in the piping. In spite of the large strides made by computational hardware 
technology, it is still extremely expensive to analyze the entire test loop utilizing high-fidelity LES 
simulations that resolve most of the scales of unsteadiness. System modeling tools on the other hand, 
are either limited in scope such as studies reported on transient pipe flows that neglect components such 
as orifices, valves and the instabilities associated with them, or oversimplify the physical mechanisms 
inherent in these components. There are two approaches used to model system components that can be 
used in conjunction with high-fidelity CFD. The first approach uses a Lumped Element Model (LEM) 
which is based on the representation of the governing physical parameters (flow coefficient, structural 
compliance) of the components by equivalent electronic analogues (resistance, capacitance) in an 
electrical circuit [7], The second approach assigns a transfer matrix (TM) to each component that relates 
inflow perturbations with outflow responses. The TM for a particular component can be derived 
mathematically from physical theory governing the working of a component, experimental or testing data, 
or computational analysis. TM techniques have been applied successfully to acoustic problems [8], and 
cavitation instability problems with turbopumps [9] and are a viable option for simulating the effects of 
instabilities in ducts, diffusers and flow control elements in rocket propulsion testing facilities. In the next 
few sections we discuss the TM approach and derive representative forms of the transfer matrix for 
common components seen in test facilities such as ducts, bends and orifices. We examine a Kelvin 
Helmholtz instability using high fidelity CFD and subsequently derive an equivalent TM from the CFD 
computations. Lastly, we discuss the propagation of an instability in a scaled down test loop which 
comprises of a duct, an orifice, a bend and a diffuser using the TM approach. 

2 

Approved for public release; distribution is unlimited. 

Copyright © 201 0 by the authors. Published by JANNAF with permission. All rights reserved. 


RELEASED - Printed documents may be obsolete; validate prior to use. 



TRANSFER MATRIX APPROACH 


A transfer matrix (TM), also known as a transmission matrix, T-matrix, or ABCD matrix, is 
traditionally employed in electrical and mechanical vibration theory [10]. In both fields, the analysis of 
complicated systems has been greatly simplified by the use of this method. Complex circuits or structures 
can be simply regarded as a “two port network” with input and output ports, as shown in Figure 1 . Matrix 
algebra can then be applied for analysis of such networks by linking the dynamic behavior between inputs 
and outputs through the coefficients of the matrix. The interconnection of parallel, series-parallel, and 
parallel-series network combinations can thus be handled by simple linear addition or multiplication of the 
transfer matrix. Since 1930, the transfer matrix method has also enjoyed considerable use in acoustic 
problems [11]. In the following sub-sections, we will develop the transfer matrices for the wave 
propagation through a uniform or variable-area duct with or without incompressible or compressible flow, 
an orifice, and a 90-degree bend, which comprises the test loop as shown in Figure 2. 
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Figure 1. A Two-Port Network. 



Consists of Ducts, an Orifice, and a 90° 
Bend. 


TRANSFER MATRIX FOR A DUCT ELEMENT 


The transfer matrix for a straight duct was derived for a wave propagating through the duct with 
and without a mean flow in a viscous medium that could either consist of an incompressible or 
compressible fluid. The pressure perturbation in the pipe is assumed to be smaller than the product of the 
density and square of the characteristic speed of sound, thereby permitting large pressure fluctuations, an 
important requirement for liquid rocket engine systems and test facilities. Furthermore, the transfer matrix 
accounts for viscous/thermal losses as well as the absorption of low frequency sound waves in the 
boundary layer of the pipe system. The detailed derivations and the theory are discussed in detail in [12]. 
The derived transfer matrix representation of the straight duct in its general form is given as 
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where k c =[fc 0 - j a+CM 0 ]/ 1 -M 0 2 where a and Z, represent the effects of viscous/thermal losses, k 0 is 

the incident wave number, I is the length of the duct and MO is the characteristic Mach Number. The 
effects of viscous friction, heat conduction and acoustic attenuation in the boundary layer included in the 
derivation here are extensively dealt in Refs. [13], [14]. Under its various approximations, the general 
form of transfer matrix can be reduced as follows: 
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COMPARING ANALYTICAL AND EXPERIMENTAL RESULTS 


(3) 

(4) 

(5) 


In this section, the comparison between the transfer matrix derivation for the 1-D wave 
propagation within a duct with mean flow and the experimental data is demonstrated. The experimental 
results are presented in [15]. For demonstration, a rectangular duct, as shown in Figure 3, is considered. 
The length of the duct is 1= 0.28 m and the width and height are 0.12 m and 0.05 m respectively. The 

duct is driven by a speaker upstream and details for the experimental setup and acoustic measurement 
are shown in [15]. The mean flow velocities are set in the experiment to be m 0 = 0, 0.1, 0.2, and 0.3. 

The comparison is made for the case of hard-walls filling the openings for liner, as shown in Figure 3. 
The comparison results for a variety of Mach Numbers at an excitation frequency of 1500 Hz are shown 
in Figure 4. The results indicate that the TM predictions agree very well with the experimental data. 

Acoustic 



Figure 3. Illustration of the Experimental Facility and the Rectangular Duct as Defined in Ref. [15]. 
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M =0.3 


Figure 4. The Acoustic Pressure at Sampling Points along the Central X-axis in a Rectangular 
Duct for various Mach numbers and f = 1500 Hz. Experimental data is compared against analytical 

TM formulations. 


TRANSFER MATRIX FOR A VARIABLE AREA DUCT ELEMENT SUCH AS DIFFUSERS AND NOZZLES 

The transfer matrix for a variable area duct/diffuser/nozzle is an extension of the one derived for 
the straight duct. However there are difficulties associated with the variation of the mean flow variables 
p 0 ,p 0 ,U 0 area change. In such a case, the problem of finding a transfer matrix becomes complex 

because the acoustic differential equations (i.e., dp' /dx and du'/dx) have variable coefficients. Two 
approaches can be considered in extending the transfer matrix for the variable area duct. In the first 
approach [16], the variable area duct was defined as an accumulation of segmented straight ducts each 
with a slightly different area (see Figure 5(a)). In the second approach also discussed in Ref. [16], the 
profile of the variable area duct/nozzle/diffuser section was defined mathematically such as with an 
exponential function and each section (See Figure 5(b)) was assumed to be comprised of constant 
coefficients of acoustic differential equations (i.e., dp'/dx and du'/dx). The segment transfer matrices 
are then computed, and these matrices are multiplied sequentially to obtain the overall transfer matrix. 


Segment Actual duct 




(a) Compilation of Constant Area Segments (b) Variable Area for each Segment 

Figure 5. Segmentation Strategies for Variable Area Duct. 
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The first step is to derive the transfer matrix for each section with area change. This can be 
written as: 
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where all entries of the matrix B are as follows 
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The next step is to divide the duct into a number of segments. The size of the segment is 
selected so that the Mach number can be assumed constant in each segment and the cross-sectional 
area profile can be locally approximated as by an exponential variation as an example, in other words 


M 0 = const. , 
and 

A = Ae™ . 

Hence, the gradient of the logarithmical variation of area is 
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Consequently, all entries of the matrix B are independent of x in each segment. The solution of 
Eqn. (6), i.e., DY/dx = BY , is given by 
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where Ax is the length of the / th segment. The value of e can be evaluated as 


? b&x =I + BAx + 


B Ax 


21 

So, the transfer matrix of the each segment is given be 
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The overall transfer matrix of wave propagation in a variable area duct is thus found from 


( 14 ) 
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( 17 ) 


r=m 

i 

For the sake of completeness, the derived matrix formulation was verified for the degenerate 
case of a straight duct and simple case for the special case of a conical pipe. 


TRANSFER MATRIX FOR AN ORIFICE 

A transfer matrix was also derived for 
an orifice-pipe system. The derivation assumed 
the mean flow to be uniform in the upstream 
duct, transition into a jet through the vena 
contracta, followed by a turbulent mixing region 
immediately downstream of the orifice, with 
flow recovery further downstream. The effect of 
a propagating perturbation is derived by 
applying the continuity and momentum 
equations in the turbulent mixing region 
downstream of the orifice, and defining the 
velocity and pressure variables as a 
combination of fluctuating and average 
quantities [17]. Subsequent linearization of the 
equations results in the coefficients of the transfer matrix. 

The transfer matrix for the orifice can be written as 

p 2 _ 1 — p 0 c 0 M 0 /3 

Ur. 0 l 


cv 



Figure 6. Schematic of a Quasi-Steady Flow 
Through an Orifice within a Duct. 


Pi 

u. 


(18) 


where the coefficient, p= S p /Sj -1 , represents the combination of the effect of the geometry and the 

vena contracta. It should be noted that the transfer matrix derived here does not account for the Kelvin 
Helmholtz instability initiated at the lip of the orifice. Extension of the TM to include the instability will be 
carried out in the future. 


TRANSFER MATRIX FOR AN TURNING DUCT 


Using a similar procedure as the orifice pipe system, a transfer matrix has also been derived for a 
disturbance/fluctuation passing through a sharp bend or a turning duct at low Strouhal Numbers [16]. The 
approach assumes that the flow separates behind the bend and the separated flow region acts as a 
blockage creating a natural vena contracta. The flow through the bend accelerates due to the reduced 
area formed by the vena contracta leading to a turbulent mixing region before the flow settles at the exit. 
Given the similarity to the flow dynamics of the orifice-duct system an analogous matrix was derived for 
sharp bends (radius of curvature to pipe diameter ratio (r/D) of less than 2) and turning ducts. The main 
distinction between the two matrices lies in the definition of the vena contracta. In the case of the orifice, 
the vena contracta is geometrically/structurally defined, whereas in the turning duct it is a function of flow 
characteristics such as Reynolds number as well as the curvature of the bend. The wall friction in the 
turbulent region is negligible 

Hence, the transmission matrix for the wave propagating through the bend is given by 
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Note that when the Mach number of the mean flow is larger than 0.3, the sharp bend is better 
described by a compressible quasi-steady flow theory [18]. The coefficient j3 depends on the extent of 
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the separation region (width), which in turn depends on the geometry of the bend and the characteristic 
Reynolds Number of the flow. 

In the case of the smooth bend ( r > 2D ), as shown in Figure 8, it is not necessary to assume a 
free jet formation with a vena contracts at the bends. Thus, it may model the smooth bend as straight 
duct [18], as discussed in the previous section on the transfer matrix for straight ducts. 


Uniform inflow 


■Pi 

■ u. 


Pi 


8 % 

P 2 h 


Jet region 
with vena 
contracta 


Turbulent 

Mixing 

region 

Uniform 

outflow 


Figure 7. Schematic of a quasi-steady flow 
through a bend within a duct. 



through a smooth (r>D) bend within a duct. 


UNSTEADY CFD ANALYSES OF THE ORIFICE WITH HYBRID RANS/LES APPROACH 


An orifice is routinely used in rocket engine test facilities to step down the pressure in a finite 
length of pipe and for flow control purposes. In this section, we discuss unsteady hybrid RANS-LES 
simulations of a pump discharge orifice in a stand-alone mode i.e. the orifice will not be coupled to any 
other flow components. The primary aim of this simulation is to characterize the instability that sets up in 
the orifice, understand the physics governing the unsteadiness and identify the unsteady modes related 
to the operation of the orifice. 

Three-dimensional simulations were carried out of a pipe with a 4-inch diameter and an orifice 
with a vena contracta of 2.08 inches. The orifice is not a sharp edge slot type orifice but rather has a 
rounded shape to guide the flow smoothly to the vena contracta. This pipe-orifice configuration has been 
used for testing at NASA SSC with liquid oxygen as the operating fluid. The upstream pipe is long enough 
to ensure that the flow is well-developed before it encounters the orifice. The flow in the pipe has a mean 
velocity of 27.64 m/s, a temperature of 91.67 K and a characteristic Reynolds Number of approximately 
4.2 million based on the pipe diameter. The time averaged solutions of pressure and velocity on the 
vertical plane are shown in Figure 9. The velocity distribution shows the flow accelerating through the 
orifice with a well defined jet extending to about a pipe diameter downstream of the orifice. Further 
downstream, the jet diffuses out to fill the pipe in the radial direction. The pressure distribution shows the 
pressure drop precipitously in the vena contracta and then recover about one diameter downstream of the 
orifice. While the drop in pressure on average is about 400 psi downstream of the orifice, a small pocket 
of very low pressure is formed in the vena contracta which may be susceptible to cavitation. 
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Figure 9. Time averaged solution of velocity and pressure distribution on the symmetry plane. 

A series of snapshots of instantaneous velocity distribution is shown in 

Figure 10 which indicates the formation of the primary jet as the flow accelerates to negotiate the 
orifice. This jet is representative of a high Reynolds Number flow due to the low viscosity of liquid 
oxygen. The snapshots of velocity distribution show the periodic generation of vortices from the lip of the 
orifice, the migration of the vortex structures into the downstream pipe section and subsequent dissipation 
of these structures further downstream in the pipe section. As a consequence of the production and 
dissipation of the vortex structures, and the dynamic interaction of the resultant jet with the wall, the 
length of the jet core is seen to periodically oscillate. The snapshots of pressure distribution are seen in 
Figure 11 which show spatially periodic patterns of low pressure (Figure 12) associated with the vorticity 
production at the lip of the orifice. The resulting Kelvin Flelmholtz type instability in the pipe/orifice system 
leads to the unsteady shedding and the periodic formation of pockets of low pressure. Pressure traces 
from numerical probes located downstream of the orifice (after the jet expands radially in the pipe) show 
two important modes that are excited: (i) a high frequency mode of 180 Hz that corresponds to the Kelvin 
Helmholtz instability and (ii) a lower frequency mode of 45 Hz with an overtone at 90 Hz. The lower 
frequency mode is an axial mode that corresponds to the interaction of the jet with the pipe wall. This 
process is illustrated in the sequence of plots in Figure 13 where first we see the jet expanding radially in 
the pipe downstream of the orifice. Since the jet expands in the confined space of the pipe, the jet 
impinges on the pipe wall, and travels upstream toward the orifice creating a feedback loop, as the 
reverse flow interacts with the primary jet leading to a sub-harmonic of the primary Kelvin Helmholtz 
instability. The lower frequency mode also manifests itself as an axial instability causing the jet core to 
fluctuate in length. 


Time 0.102208 


Time 0.132208 



Time 0.1 26208 
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Figure 10. Snapshots of Instantaneous Axial Velocity Distribution. 
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Figure 11. Instantaneous Snapshots of Pressure Distribution on the Vertical Plane. 
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Figure 12. Pressure Spectra Showing 
Dominant Modes. 



Interaction with wall creates a feedback loop which in turn couples 
with a sub-harmonic of the primary instability 



The secondary lower frequency instability also manifests 
itself as a axial instability 


Figure 13. Snapshots of Axial Velocity Distribution 
Showing the Evolution of the Instability Modes in the 
Computational Domain. 


COMPUTATION OF TRANSFER MATRIX FROM CFD SIMULATIONS 


Here we discuss the procedure used to compute the transfer matrix of a particular component 
directly from the CFD simulations. In the case of the orifice simulations, the history of pressure and 
velocity fluctuation was recorded by numerical probes at a series of points upstream and downstream of 
the orifice. For the purpose of demonstration we consider the fluctuation history at a point one-half inch 
upstream of the orifice and nine inches downstream of the orifice (See Figure 14). It is important to note 
here that (i) the downstream point is in the turbulent mixing region after the flow has re-attached to the 
wall and (ii) the fluctuations that are used in the derivation of the transfer matrix should be averaged 
across the cross-section plane rather than use isolated points as has been done in this demonstration 
case. The pressure and velocity spectra are shown in Figure 15 for both the points on the upstream and 
downstream sides. 
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x = 9 


x= -0.5 

Figure 14. Orifice Domain Showing location of the Upstream Probe(31) at x=-0.5 
and Downstream Probe(1 5) located at x=9. 


Upstream Side Downstream Side 



Frequency (Hz) Frequency (Hz) 

Figure 15. Pressure and velocity spectra at the upstream and downstream probe points. 


The amplitudes of the pressure and velocity at the low frequency (46 Hz) and high frequency (180 
Hz) for both points are extracted from their respective spectra and used to determine the coefficients of 
the transfer matrix given by the system below. As was shown previously in the discussion of the transfer 
matrix the subscript 2 stands for the fluctuations at the downstream end and the subscript 1 refers to the 
fluctuations at the upstream end. The amplitudes at the two dominant frequencies from the CFD 
simulations are used to solve the linear system of equations, thereby obtaining the coefficients of the 
transfer matrix. Both the pressure and velocity fluctuations downstream of the orifice show relatively 
greater dependence on the velocity fluctuation upstream of the orifice. This is also seen in the 
analytically derived expressions for the orifice that were discussed previously. 


where 
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SIMULATION OF A TEST LOOP WITH MULTIPLE COMPONENTS 


In this task we demonstrate the ability of the instability analysis framework to analyze a scaled 
down test loop with more than one component. The scaled down facility is shown in Figure 16, and 
comprises of a duct component, an orifice, a 90° bend, and a conical expansion (diffuser). The pipe 

system has input pressure and streamwise velocity fluctuation magnitudes \p[\ = 200 psi and 
\u[\ = 5 m/s , respectively. Other parameters are given for the liquid fluid as p 0 — lOOOkg/m 3 , 
c 0 = 1200 ml s , mean input pressure P 0 = 800 psi and mean input velocity U 0 =20 m/s . 



T 

4' 

1 


Pa 


Ps 



Figure 16 . A schematic of a scaled down system with multiple components (drawing not to scale). 

In order to derive the TM representation of the demo pipe system, it is assumed that the orifice is 
located in the middle of the orifice pipe section, as shown in Figure 16. Furthermore, the orifice 
component is assumed to occupy two pipe diameters upstream (i.e., Duct B) and two pipe diameters 
downstream (Duct C). The length of the downstream Duct C, in particular, depends upon the length of 
the mixing region. For the bend section, because it is a smooth bend ( r > 3D), it is thus treated as a 
duct component as explained earlier. The physical motivation for this comes from the lack of a significant 
separation bubble size for smooth bends but, again, is expected to be dependent on flow Reynolds 
number and surface roughness. This model can be fine tuned with the aid of RANS CFD simulations 
and/or experiments. Consequently, the TM representation for this scaled down system is given by a 
serial product of the respective transfer matrices of each component 
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Figure 17 shows the TM prediction for the amplitude of the perturbation pressure and velocity at 
the outlet of the duct A p/,u / and the expansion p 5 ' ,u 5 ' versus frequency, as indicated in Figure 

16. In this case, the pressure p[ and velocity u[ perturbations are in phase. The amplitude of the 
perturbation pressure and velocity are normalized by | p[\ and |m[| , respectively. The results indicate that 

at certain frequencies, the perturbation pressure within the system at a given location can reach a local 
maximum. In other words, the system (section) is resonant at these frequencies, which could lead to 
large fluctuating forces on the piping system. For instance, Duct A is resonant at 100 Hz for the given 
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conditions. Figure 18 shows the TM prediction for a similar case when p[ and u[ are out of phase by 

90°. The change of the excitation condition, which is dependent on components upstream of Duct A, 
shifts slightly the resonant frequency(ies) of the system (section). 

This model has multiple benefits. These results can be used to predict potential frequencies of 
oscillation or resonance. If the resonant frequencies are near the estimated resonant frequencies of the 
piping system (obtained from a modal structural analysis), then the geometry of system components can 
be adjusted to “detune” the respective frequencies away from each other. Furthermore, as additional 
data (experimental or CFD) are compiled, the TM model can be adjusted to improve its accuracy. 



Figure 17. Perturbation pressure and velocity 
within a demo pipe system with p[ and u[ in 
phase. 



Figure 18. Perturbation pressure and velocity 
within a demo pipe system with p[ and u[ out 
of phase by 90°. 


SUMMARY AND CONCLUSIONS 

In summary, the significance of this paper is the development of a multi-level tool that can resolve 
a governing instability in a liquid rocket propulsion system or test facility, and identify modes that amplify 
as the instability travels through the various components that comprise the system. The primary source 
of the unsteadiness is modeled with a high-fidelity hybrid RANS/LES based CFD methodology that has 
been previously used to study instabilities in feed systems, such as in orifices/venturies and transient 
valve operation in test facilities. System response to the instability can be predicted through the use of a 
transfer matrix approach where the individual responses of all the components in the system can be 
integrated together in real-time to identify potentially harmful resonant modes. The resultant multi-level 
tool can be used to mitigate these resonant modes by providing design support for (i) alternate system 
configurations and/or (ii) control mechanisms for suppression of the resonant modes. 

The focus of this paper relates to the development of a system level tool based on the Transfer 
Matrix approach for the evaluation of a propagating instability. The transfer matrix is used to link the 
pressure and velocity fluctuations at the inlet of a given component with the analogous fluctuations at the 
exit of the component. The coefficients of the transfer matrix are derived rigorously from one-dimensional 
wave theory principles based on the dominant flow physics characteristic of the component. Transfer 
Matrices were derived for components that are commonly seen in rocket engine systems/test stands such 
as pipe systems that consist of straight ducts and 90° bends, flow control elements, such as orifices, and 
variable area ducts, such as nozzles and diffusers. One of the major advantages of the transfer matrix 
approach is that the propagation of a disturbance through a system loop can be analyzed by assembling 
a simple product of the transfer matrices of the components that comprise the system. An integral part of 
this paper involved the demonstration of the propagation of an instability through a system that consisted 
of a straight duct, an orifice, a turning duct and an exit diffuser. 
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One of the important goals of the program was to utilize high-fidelity CFD to understand the 
important flow physics responsible for the governing instability and to characterize the important modes 
associated with hydrodynamic instability. Having characterized the instability with high-fidelity CFD, the 
next goal was to determine a procedure for utilizing this information with the transfer matrix approach 
used for propagating the instability through all the components that comprise a propulsion system loop. 
To this end, three-dimensional simulations with the hybrid RANS-LES model were carried out for an 
orifice operating with liquid oxygen. Simulations of the orifice show periodic shedding from the lip of the 
orifice resulting from the onset of the Kelvin Helmholtz instability and the important modes were identified 
as part of the simulation. Furthermore, information from the CFD simulations was used to construct the 
coefficients of a transfer matrix. This transfer matrix can be embedded in the system level tool discussed 
in the paper to investigate dynamic behavior of rocket engine systems and test facilities. 
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